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Abstract Phase oscillators are a common starting point for the reduced description 
of many single neuron models that exhibit a strongly attracting limit cycle. The frame- 
work for analysing such models in response to weak perturbations is now particularly 
well advanced, and has allowed for the development of a theory of weakly connected 
neural networks. However, the strong-attraction assumption may well not be the nat- 
ural one for many neural oscillator models. For example, the popular conductance 
based Morris-Lecar model is known to respond to periodic pulsatile stimulation in a 
chaotic fashion that cannot be adequately described with a phase reduction. In this pa- 
per, we generalise the phase description that allows one to track the evolution of dis- 
tance from the cycle as well as phase on cycle. We use a classical technique from the 
theory of ordinary differential equations that makes use of a moving coordinate sys- 
tem to analyse periodic orbits. The subsequent phase -amplitude description is shown 
to be very well suited to understanding the response of the oscillator to external stim- 
uli (which are not necessarily weak). We consider a number of examples of neural 
oscillator models, ranging from planar through to high dimensional models, to illus- 
trate the effectiveness of this approach in providing an improvement over the standard 
phase-reduction technique. As an explicit application of this phase-amplitude frame- 
work, we consider in some detail the response of a generic planar model where the 
strong-attraction assumption does not hold, and examine the response of the system 
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to periodic pulsatile forcing. In addition, we explore how the presence of dynamical 
shear can lead to a chaotic response. 

Keywords Phase-amplitude • Oscillator • Chaos • Non-weak coupling 

List of Abbreviations 

ML Morris-Lecar 
FHN FitzHugh-Nagumo 
CS Connor-Stevens 
LE Lyapunov exponent 

1 Introduction 

One only has to look at the plethora of papers and books on the topic of phase oscil- 
lators in mathematical neuroscience to see the enormous impact that this tool from 
dynamical systems theory has had on the way we think about describing neurons and 
neural networks. Much of this work has its roots in the theory of ordinary differen- 
tial equations (ODEs) and has been promoted for many years in the work of Winfree 
[1], Guckenheimer [2], Holmes [3], Kopell [4], Ermentrout [5] and Izhikevich [6] to 
name but a few. For a recent survey, we refer the reader to the book by Ermentrout 
and Terman [7]. At heart, the classic phase reduction approach recognises that if a 
high dimensional non-linear dynamical system (as a model of a neuron) exhibits a 
stable limit cycle attractor then trajectories near the cycle can be projected onto the 
cycle. 

A natural phase variable is simply the time along the cycle (from some arbitrary 
origin) relative to the period of oscillation. The notion of phase can even be extended 
off the cycle using the concept of isochrons [1]. They provide global information 
about the 'latent phase', namely the phase that will be asymptotically returned to for 
a trajectory with initial data within the basin of attraction of an exponentially stable 
periodic orbit. More technically, isochrons can be viewed as the leaves of the invariant 
foliation of the stable manifold of a periodic orbit [8]. In rotating frame coordinates 
given by phase and the leaf of the isochron foliation, the system has a skew-product 
structure, i.e. the equation of the phase decouples. However, it is a major challenge 
to find the isochron foliation, and since it relies on the knowledge of the limit cycle it 
can only be found in special cases or numerically. There are now a number of com- 
plementary techniques that tackle this latter challenge, and in particular we refer the 
reader to work of Guillamon and Huguet [9] (using Lie symmetries) and Osinga and 
Moehlis [10] (exploiting numerical continuation). More recent work by Mauroy and 
Mezic [11] is especially appealing as it uses a simple forward integration algorithm, 
as illustrated in Fig. 1 for a Stuart-Landau oscillator. However, it is more common 
to side-step the need for constructing global isochrons by restricting attention to a 
small neighbourhood of the limit cycle, where dynamics can simply be recast in the 
reduced form 0 = 1, where 0 is the phase around a cycle. This reduction to a phase 
description gives a nice simple dynamical system, albeit one that cannot describe evo- 
lution of trajectories in phase- space that are far away from the limit cycle. However, 
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Fig. 1 Isochrons of a 
Stuart-Landau oscillator model: 
i = Xx/2 — (Xc/2 + co)y — 
X(x 2 + y 2 )(x-cy)/2, 
y = (Xc/2 + co)x+Xy/2- 
X(x 2 + y 2 )(cx + y)/2. The 
black curve represents the 
periodic orbit of the system, 
which is simply the unit circle 
for this model. The blue curves 
are the isochrons obtained using 
the (forward) approach of 
Mauroy and Mezic [11]. The 
green dots are analytically 
obtained isochronal points [15]. 
Parameter values are X = 2.0, 
c = 1 .0, and co =1.0 

the phase reduction formalism is useful in quantifying how a system (on or close 
to a cycle) responds to weak forcing, via the construction of the infinitesimal phase 
response curve (iPRC). For a given high dimensional conductance based model this 
can be solved for numerically, though for some normal form descriptions closed form 
solutions are also known [12]. The iPRC at a point on cycle is equal to the gradient 
of the (isochronal) phase at that point. This approach forms the basis for construct- 
ing models of weakly interacting oscillators, where the external forcing is pictured 
as a function of the phase of a firing neuron. This has led to a great deal of work 
on phase-locking and central pattern generation in neural circuitry (see, for example 
[13]). Note that the work in [9] goes beyond the notion of iPRC and introduces in- 
finitesimal phase response surfaces (allowing evaluation of phase advancement even 
when the stimulus is off cycle), and see also the work in [14] on non-infinitesimal 
PRCs. 

The assumption that phase alone is enough to capture the essentials of neural re- 
sponse is one made more for mathematical convenience than being physiologically 
motivated. Indeed, for the popular type I Morris-Lecar (ML) firing model with stan- 
dard parameters, direct numerical simulations with pulsatile forcing show responses 
that cannot be explained solely with a phase model [16]. The failure of a phase de- 
scription is in itself no surprise and underlies why the community emphasises the 
use of the word weakly in the phrase "weakly connected neural networks". Indeed, 
there are a number of potential pitfalls when applying phase reduction techniques 
to a system that is not in a weakly forced regime. The typical construction of the 
phase response curve uses only linear information about the isochrons and non-linear 
effects will come into play the further we move away from the limit cycle. This prob- 
lem can be diminished by taking higher order approximations to the isochrons and 
using this information in the construction of a higher order PRC [17]. Even using 
perfect information about isochrons, the phase reduction still assumes persistence of 
the limit-cycle and instantaneous relaxation back to cycle. However, the presence of 
nearby invariant phase- space structures such as (unstable) fixed points and invari- 
ant manifolds may result in trajectories spending long periods of time away from 
the limit cycle. Moreover, strong forcing will necessarily take one away from the 
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neighbourhood of a cycle where a phase description is expected to hold. Thus, de- 
veloping a reduced description, which captures some notion of distance from cycle 
is a key component of any theory of forced limit cycle oscillators. The development 
of phase-amplitude models that better characterise the response of popular high di- 
mensional single neuron models is precisely the topic of this paper. Given that it is 
a major challenge to construct an isochronal foliation we use non-isochronal phase- 
amplitude coordinates as a practical method for obtaining a more accurate description 
of neural systems. Recently, Medvedev [18] has used this approach to understand in 
more detail the synchronisation of linearly coupled stochastic limit cycle oscillators. 

In Sect. 2, we consider a general coordinate transformation, which recasts the dy- 
namics of a system in terms of phase- amplitude coordinates. This approach is directly 
taken from the classical theory for analysing periodic orbits of ODEs, originally con- 
sidered for planar systems in [19], and for general systems in [20]. We advocate it 
here as one way to move beyond a purely phase-centric perspective. We illustrate the 
transformation by applying it to a range of popular neuron models. In Sect. 3, we 
consider how inputs to the neuron are transformed under these coordinate transfor- 
mations and derive the evolution equations for the forced phase- amplitude system. 
This reduces to the standard phase description in the appropriate limit. Importantly, 
we show that the behaviour of the phase- amplitude system is much more able to cap- 
ture that of the original single neuron model from which it is derived. Focusing on 
pulsatile forcing, we explore the conditions for neural oscillator models to exhibit 
shear induced chaos [16]. Finally in Sect. 4, we discuss the relevance of this work to 
developing a theory of network dynamics that can improve upon the standard weak 
coupling approach. 

2 Phase-Amplitude Coordinates 

Throughout this paper, we study the dynamics prescribed by the system x = f(x), 
x G R", with solutions x = x(t) that satisfy x(0) = xo eW 1 . We will assume that 
the system admits an attracting hyperbolic periodic orbit (namely one zero Floquet 
exponent and the others having negative real part), with period A, such that u(t) = 
u(t + A). A phase 6 e [0, A) is naturally defined from 6(u(t)) = £modA. It has 
long been known in the dynamical systems community how to construct a coordinate 
system based on this notion of phase as well as a distance from cycle; see [20] for 
a discussion. In fact, Ermentrout and Kopell [21] made good use of this approach to 
derive the phase-interaction function for networks of weakly connected limit-cycle 
oscillators in the limit of infinitely fast attraction to cycle. However, this assumption 
is particularly extreme and unlikely to hold for a broad class of single neuron models. 
Thus, it is interesting to return to the full phase-amplitude description. In essence, the 
transformation to these coordinates involves setting up a moving orthonormal system 
around the limit cycle. One axis of this system is chosen to be in the direction of the 
tangent vector along the orbit, and the remaining are chosen to be orthogonal. We 
introduce the normalised tangent vector § as 
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Fig. 2 Demonstration of the 
moving orthonormal coordinate 
system along an orbit segment. 
As t evolves from t$ to t\ , the 
coordinates vary smoothly. In 
this planar example, £ always 
points to the outside of the orbit 



^(6(t 0 )) 

^O(t 0 )) 




The remaining coordinate axes are conveniently grouped together as the columns of 
an n x (n — 1) matrix f . In this case we can write an arbitrary point x as 



x(0,p) = u(0) + S(0)p. 



(2) 



Here, |p| represents the Euclidean distance from the limit cycle. A caricature, in R 2 , 
of the coordinate system along an orbit segment is shown in Fig. 2. Through the use 
of the variable p, we can consider points away from the periodic orbit. Rather than 
being isochronal, lines of constant 0 are simply straight lines that emanating from a 
point on the orbit in the direction of the normal. The technical details of specifying 
the orthonormal coordinates forming f are discussed in Appendix A. 

Upon projecting the dynamics onto the moving orthonormal system, we obtain the 
dynamics of the transformed system: 



= l + /l(0,p), 



p = A(6)p + f 2 (0, p), 



where 



fl(P, P ) = 

f2(0,p) = 

h(0,p) = 



-h T (0, p)^p + h T (0, p)[f(u + fp) - /(«)], 
~^% pfl + ^[ f(U + " f{u) ~ D/ ^' 



du 
dO 



+ ? d0 P 



A(0) = t; T 



d(9 JS 



(3) 

(4) 
(5) 

(6) 



and D/ is the Jacobian of the vector field / evaluated along the periodic orbit u. 
The derivation of this system may be found in Appendix A. It is straightforward to 
show that /i(0, p) -> 0 as |p| -> 0, / 2 (0, 0) = 0 and that 3/ 2 (6>, 0)/3p = 0. In the 
above, f\ captures the shear present in the system, that is, whether the speed of 6 
increases or decreases dependent on the distance from cycle. A precise definition for 
shear is given in [22]. Additionally, A (6) describes the 0 -dependent rate of attraction 
or repulsion from cycle. 

It is pertinent to consider where this coordinate transformation breaks down, 
that is, where the determinant of the Jacobian of the transformation K = det[dx/d0 
dx/dp] vanishes. This never vanishes on-cycle (where p = 0), but may do so for 
some |p | = k > 0. This sets an upper bound on how far away from the limit cycle we 
can describe the system using these phase- amplitude coordinates. In Fig. 3, we plot 
the curve along which the transformation breaks down for the ML model. We observe 
that, for some values of 0 , k is relatively smaller. The breakdown occurs where lines 
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Fig. 3 This figure shows the 
determinant K of the 
phase-amplitude transformation 
for the ML model. Colours 
indicate the value of K. The red 
contour indicates where K = 0, 
and thus where the coordinate 
transformation breaks down. 
Note how all of the values for 
which this occurs have p < 0. 
Parameter values as in 
Appendix B.l 
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of constant 6 cross, and thus the transformation ceases to be invertible, and these val- 
ues of 6 correspond to points along which the orbit has high curvature. We note that 
this issue is less problematical in higher dimensional models. 
If we now consider the driven system, 

x = f(x) + sg(x,t), xeR n , (7) 

where s is not necessarily small, we may apply the same transformation as above to 
obtain the dynamics in (0, p) coordinates, where 6 e [0, A) and p e W 1-1 , as 

0 = 1 + /i(0, p) + sh T (0, p)g(u(0) + ?(0)p, t), (8) 
p = A(0)p + f 2 (0, p) + st; T B(0, p)g{u(0) + ?(0)p, t), (9) 

with 

B(0,p)=I n -^ph T (0,p), (10) 
aO 

and l n is the n x n identity matrix. Here, h and B describe the effect in terms of 0 
and p that the perturbations have. Details of the derivation are given in Appendix A. 
For planar models, B = I2. To demonstrate the application of the above coordinate 
transformation, we now consider some popular single neuron models. 



2.1 A 2D Conductance Based Model 



The ML model was originally developed to describe the voltage dynamics of barnacle 
giant muscle fibre [23], and is now a popular modelling choice in computational 
neuroscience [7]. It is written as a pair of coupled non-linear ODEs of the form 

Cv = I{t) - g\(v - V\) - gKW(v - Vk) ~ gCa^ocO)0 V Csl ), 

(11) 

w = 0(woc(u) - w)/t w (v). 

Here, v is the membrane voltage, whilst w is a gating variable, describing the frac- 
tion of membrane ion channels that are open at any time. The first equation expresses 
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Fig. 4 Typical trajectory of the ML model of the transformed system. Left: Time evolution of 6 and p. 
Right: Trajectory plotted in the (v, w) phase plane. We see that when p has a local maximum, the evolution 
of 6 slows down, corresponding to where trajectories pass near to the saddle node. Parameter values as in 
Appendix B.l 



Kirchoff 's current law across the cell membrane, with I(t) representing a stimulus in 
the form of an injected current. The detailed form of the model is completed in Ap- 
pendix B.l. The ML model has a very rich bifurcation structure. Roughly speaking, 
by varying a constant current I(t) = Io, one observes, in different parameter regions, 
dynamical regimes corresponding to sinks, limit cycles, and Hopf, saddle-node and 
homoclinic bifurcations, as well as combinations of the above. These scenarios are 
discussed in detail in [7] and [24] . 

As the ML model is planar, p is a scalar, as are the functions A and f\^. This 
allows us to use the moving coordinate system to clearly visualise parts of phase 
space where trajectories are attracted towards the limit cycle, and parts in which they 
move away from it, as illustrated in Fig. 4. The functions and A, evaluated at 
p = —0.1 are shown in Fig. 5. The evolution of 0 is mostly constant, however we 
clearly observe portions of the trajectories where this is slowed, along which p ~ 0. 
In fact, this corresponds to where trajectories pass near to the saddle node, and the 
dynamics stall. This occurs around 0 = 17, and in Fig. 5 we see that both A(0) and 
f\(0, p) are indeed close to 0. The reduced velocities of trajectories here highlights 
the importance of considering other phase space structures in forced systems, the 
details of which are missed in standard phase only models. Forcing in the presence 
of such structures may give rise to complex and even chaotic behaviours, as we shall 
see in Sect. 3. 

In the next example, we show how the same ideas go across to higher dimensional 
models. 

2.2 A 4D Conductance Based Model 

The Connor-Stevens (CS) model [25] is built upon the Hodgkin-Huxley formalism 
and comprises a fast Na + current, a delayed K + current, a leak current and a transient 
K + current, termed the A-current. The full CS model consists of 6 equations: the 
membrane potential, the original Hodgkin-Huxley gating variables, and an activating 
and inactivating gating variable for the A-current. Using the method of equivalent 



4y Springer 



Page 8 of 22 



K.C.A. Wedgwood et al. 




potentials [26], we may reduce the dimensionality of the system, to include only 4 
variables. The reduced system is 

Cv = —F(v,u,a,b) + /, u = G(v,u), 

. X OQ (v)-X (12) 
X = , Xe{a,b], 

where 

F(v, u, a, b) = gKn^Miv - v K ) + ( gNa^oo(w)m^ ) (i;)(i; - v Na ) 

+ gil(v - vi) + g a a 3 b(v - u a ). (13) 

The details of the reduced CS model are completed in Appendix B.2. The solutions 
to the reduced CS model under the coordinate transformation may be seen in Fig. 6, 
whilst, in Fig. 7, we show how this solution looks in the original coordinates. As for 
the ML model, 0 evolves approximately constantly throughout, though this evolution 
is sped up close to 0 = A. The trajectories of the vector p are more complicated, 
but note that there is regularity in the pattern exhibited, and that this occurs with 
approximately the same period as the underlying limit cycle. The damping of the 
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Fig. 6 Solution of the transformed CS system. Top: Time evolution of 6. Bottom: Time evolution of p 
coordinates. Upon transforming these coordinates back to the original ones, we arrive at Fig. 7. Parameter 
values given in Appendix B.2. In this parameter regime, the model exhibits type I firing dynamics 

amplitude of oscillations in p over successive periods represents the overall attraction 
to the limit cycle, whilst the regular behaviour of p represents the specific relaxation 
to cycle as shown in Fig. 7. Additional file 1 shows a movie of the trajectory in 
(v,u,b) coordinates with the moving orthonormal system superimposed, as well as 
the solution in p for comparison. 

3 Pulsatile Forcing of Phase- Amplitude Oscillators 

We now consider a system with time-dependent forcing, given by (7) with 

g(x,t) = J](5a-nr),0,...,0) r , (14) 

neZ 

where 8 is the Dirac 8 -function. This describes T -periodic kicks to the voltage vari- 
able. Even such a simple forcing paradigm can give rise to rich dynamics [16]. For 
the periodically kicked ML model, shear forces can lead to chaotic dynamics as folds 
and horseshoes accumulate under the forcing. This means that the response of the 
neuron is extremely sensitive to the initial phase when the kicks occur. In terms of 
neural response, this means that the neuron is unreliable [27]. 

The behaviour of oscillators under such periodic pulsatile forcing is the subject of 
a number of studies; see, e.g. [27-30]. Of particular relevance here is [27], in which 
a qualitative reasoning of the mechanisms that bring about shear in such models is 



Springer 



Page 10 of 22 



K.C.A. Wedgwood et al. 




Fig. 7 Transformed trajectory in (v,u,b) space of the phase-amplitude description of the reduced CS 
model. The dotted black line is the phase amplitude solution transformed in the original coordinates, whilst 
the coloured orbit is the underlying periodic orbit, where the colour corresponds to the phase along the 
orbit. The solution of the phase-amplitude description of the model in (0, p) coordinates is shown in Fig. 6 



supplemented by direct numerical simulations to detect the presence of chaotic solu- 
tions. For the ML model in a parameter region close to the homoclinic regime, kicks 
can cause trajectories to pass near the saddle-node, and folds may occur as a result 
[16]. 

We here would like to compare full planar neural models to the simple model, 
studied in [27]: 

0 = l+crp, p = -kp + eP(0)^28(t-nT). (15) 

neZ 

This system exhibits dynamical shear, which under certain conditions, can lead to 
chaotic dynamics. The shear parameter a dictates how much trajectories are 'sped 
up' or 'slowed down' dependent on their distance from the limit cycle, whilst A is 
the rate of attraction back to the limit cycle, which is independent of 6 . Supposing 
that the function P is smooth but non-constant, trajectories will be taken a variable 
distance from the cycle upon the application of the kick. When kicks are repeated, 
this geometric mechanism can lead to repeated stretching and folding of phase space. 
It is clear that the larger a is in (15), the more shear is present, and the more likely 
we are to observe the folding effect. In a similar way, smaller values of A mean that 
the shear has longer to act upon trajectories and again result in a greater likelihood of 
chaos. Finally, to observe chaotic response, we must ensure that the shear forces have 
sufficient time to act, meaning that T, the time between kicks must not be too small. 
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Fig. 8 Stretch-and-fold action of a kick followed by relaxation in the presence of shear. The thin black 
lines are the isochrons of the system, which in the case of the linear model (15), are simply straight lines 
with slope —X/cr. The thin grey line at p = 0 represents the limit cycle, which is kicked, at t = to by 
P(0) = sin(#) with strength s = 1 to the solid curve. After this, the orbits are allowed to evolve under the 
flow generated by the continuous part of the system. The dashed and dotted curves represent the image 
of the kicked solid curve under this flow, at times t\ and t2, respectively. The green marker shows how 
one point, x(t$) evolves under the flow, first to x(t\) and then to xfo), following the isochron as it relaxes 
back to the limit cycle. The effect of the shear forces and the subsequent folding, caricatured by the blue 
arrows can clearly be seen 

This stretching and folding action can clearly lead to the formation of Smale horse- 
shoes, which are well known to lead to a type of chaotic behaviour. However, horse- 
shoes may co-exist with sinks, meaning the resulting chaotic dynamics would be 
transient. Wang and Young proved that under appropriate conditions, there is a set of 
T of positive Lebesgue measure for which the system experiences a stronger form of 
sustained, chaotic behaviour, characterised by the existence of a positive Lyapunov 
exponent for almost all initial conditions and the existence of a 'strange attractor'; 
see, e.g. [28-30]. 

By comparing with the phase-amplitude dynamics described by Eqs. (8)-(9), we 
see that the model of shear considered in (15) is a proxy for a more general system, 
with fi(0, p) -> ap, A(0) -> —A and h(0, p) -> 0, and f (0) -> P(0). 

To gain a deeper insight into the phenomenon of shear induced chaos, it is pertinent 
to study the isochrons of the limit cycle for the linear model (15), where the isochrons 
are simply lines with slope —A/a. In Fig. 8, we depict the isochrons and stretch and 
fold action of shear. The bold thin grey line at p = 0 is kicked, at t = to, to the bold 
solid curve, where P(0) = sin(#), as studied in [16] and this curve is allowed to 
evolve under the dynamics with no further kicks through the dashed curve at t = t\ 
and ultimately to the dotted curve at t = t2, which may be considered as evolutions of 
the solid curve by integer multiples of the natural period of the oscillator. Every point 
of the dotted curve traverses the isochron it is on at to until ^- The green marker shows 
an example of one such point evolving along its associated isochron. The folding 
effect of this is clear in the figure, and further illustrated in the video in Additional 
file 2. 
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Fig. 9 Shear induced folding in the ML model, with parameters as in Fig. 5. The red curve in all panels 
represents the limit cycle of the unperturbed system, whilst the green dotted line represents the stable 
manifold of the saddle node, indicated by the orange marker. We begin distributing points along the limit 
cycle and then apply an instantaneous kick taking v \-> v + s, where s = —2.0, leaving w unchanged. This 
essentially moves all phase points to the left, to the blue curve. The successive panels show the image 
of this set of points after letting points evolve freely under the system defined by the ML equations, and 
then apply the kick again. The curves shown are the images of the initial phase points just after each 
kick, as indicated in the figure. We can clearly observe the shear induced folding. Parameter values as in 
Appendix B.l 



This simple model with a harmonic form for P{6) provides insight into how 
strange attractors can be formed. Kicks along the isochrons or ones that map 
isochrons to one another will not produce strange attractors, but merely phase-shifts. 
What causes the stretching and folding is the variation in how far points are moved as 
measured in the direction transverse to the isochrons. For the linear system (15) vari- 
ation in this sense is generated by any non-constant P{6)\ the larger the ratio crs/X, 
the larger the variation (see [16] for a recent discussion). 

The formation of chaos in the ML model is shown in Fig. 9. Here, we plot the 
response to periodic pulsatile forcing, given by (14), in the (v, w) coordinate system. 
This clearly illustrates a folding of phase space around the limit-cycle, and is further 
portrayed in the video in Additional file 3. We now show how this can be understood 
using phase-amplitude coordinates. 

Compared to the phenomenological system (15), models written in phase- 
amplitude coordinates as (8)-(9) have two main differences. The intrinsic dynamics 
(without kicks) are non-linear, and the kick terms appear in both equations for 0 and 
p (not just p). Simulations of (8)-(9) for both the FHN and ML models, with s = 0. 1, 
show that the replacement of f\(0, p) by op, dropping f2(0, p) (which is quadratic 
in p), and setting A(0) = —X does not lead to any significant qualitative change in 
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behaviour (for a wide range of o, X > 0). We therefore conclude that, at least when 
the kick amplitude s is not too large, it is more important to focus on the form of 
the forcing in the phase- amplitude coordinates. In what follows, we are interested in 
discovering the effects of different functional forms of the forcing term multiplying 
the 8 -function, keeping other factors fixed. As examples, we choose those forcing 
terms given by transforming the FHN and the ML models into phase-amplitude co- 
ordinates. To find these functions, we first find the attracting limit cycle solution in 
the ML model (11) and FHN model (52) using a periodic boundary value problem 
solver and set up the orthonormal coordinate system around this limit cycle. Once the 
coordinate system is established, we evaluate the functions h(0, p) and B(0, p) (that 
appear in Eqs. (8) and (9)). For planar systems, we have simply that B(0, p) =h- 
Using the forcing term (14), we are only considering perturbations to the voltage 
component of our system and thus only the first component of h, and the first column 
of B will make a non-trivial contribution to the dynamics. We define P\ as the first 
component of h and P 2 as the first component of £ . We wish to force each system 
at the same ratio of the natural frequency of the underlying periodic orbit. To ease 
comparison between the system with the ML forcing terms and the FHN forcing 
terms, we rescale 0 \-> 6/ A so that 0 e [0, 1) in what follows. Implementing the 
above choices leads to 

0 = 1 + op + sPi(0, p)^8(t - nT), 

(16) 



p = -Xp + eP 2 (0) ^8{t-nT). 



It is important to emphasise that P\ 2 are determined by the underlying single neuron 
model (unlike in the toy model (15)). As emphasised in [31], one must take care 
in the treatment of the state dependent 'jumps' caused by the 8 -functions in (16) 
to accommodate the discontinuous nature of 0 and p at the time of the kick. To 
solve (16), we approximate 8(t) with a normalised square pulse 8 T (t) of the form 



«r(0 = 



0, t < 0, 

1/T, 0<t<T, (17) 

0, t > r, 



where t « 1. This means that for (n — l)T + r < t < nT, n e Z, the dynamics 
are governed by the linear system (0, p) = (1 + op, —Xp). This can be integrated 
to obtain the state of the system just before the arrival of the next kick, (0 n , p n ) = 
(0(nT), p(nT)), in the form 



0(0 + nT - t + -p(0(l - Q-^ nT -^) 
X 

-k(nT-t) 



modi, (18) 



p„ = p(0e- A ^- r \ (19) 

In the interval nT <t < nT + t and using (17) we now need to solve the system of 
non-linear ODEs 

0 = 1 + op + -Pi(0, p), p = -Xp + -P 2 (0). (20) 

T T 
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Rescaling time as t = nT + xs, with 0 < s < 1, and writing the solution (0, p) 
as a regular perturbation expansion in powers of x as (0(s), p(s)) = (9o(s) + 
r^i(5), poCO + tpifa)) + • • • , we find after collecting terms of leading order in x 
that the pair (Oo(s), po(s)) is governed by 

^ = ePi(Oo(s), Po(s)), ^ = eft(ft)0)), 0 < s < 1, (21) 

with initial conditions (#o(0), Po(0)) = (0 n ,p n ). The solution (0o(l)» Po(l)) = 
(#+,p+) (obtained numerically) can then be taken as initial data (0(t),p(t)) = 
(0+, p+) in (18)— (19) to form the stroboscopic map 



^+1 



9+ + r + ^p+(l-e-^) 



modi, (22) 

A 

p, +1 =p+e- Ar . (23) 



Note that this has been constructed using a matched asymptotic expansion, using 
(17), and is valid in the limit x —> 0. For weak forcing, where s <gi 1, Pi, 2 var Y slowly 
through the kick and can be approximated by their values at (0„ , p n ) so that to 0 (e 2 ) 



7/1+1 = 



0 n + T + £ Pi (0 n , p n ) + ^(Pn+ eP 2 (0 n ))(l - e~ kT ) 



modi, (24) 



P^+i = (pn+^2(^))e~ Ar . (25) 

Although this explicit map is convenient for numerical simulations, we prefer to work 
with the full stroboscopic map (22)-(23), which is particularly useful for comparing 
and contrasting the behaviour of different planar single neuron models with arbitrary 
kick strength. As an indication of the presence of chaos in the dynamics resulting 
from this system, we evaluate the largest Lyapunov exponent of the map (22)-(23) 
by numerically evolving a tangent vector and computing its rate of growth (or con- 
traction); see e.g. [32] for details. 

In Fig. 10, we compare the functions for both the FHN and the ML models. 
We note that P2 for the FHN model is near 0 for a large set of 6, whilst the same is 
true for Pi for the ML model. This means that kicks in the FHN model will tend to 
primarily cause phase shifts, whilst the same kicks in the ML model will primarily 
cause shifts in amplitude. 

We plot in the top row of Fig. 1 1 the pair (0 n , 0 n +\ ) , for (24)-(25) for the FHN and 
ML models. For the FHN model, we find that the system has Lyapunov exponent of 
—0.0515 < 0. For the ML model the Lyapunov exponent is 0.6738 > 0. This implies 
that differences in the functional forms of Pi ,2 can help to explain the generation of 
chaos. 

Now that we know the relative contribution of kicks in v to kicks in (0, p), it is 
also useful to know where kicks actually occur in terms of 0 as this will determine the 
contribution of a train of kicks to the (0, p) dynamics. In Figs. 11c and d, we plot the 
distribution of kicks as a function of 6 . For the ML model, we observe that the kicks 
are distributed over all phases, while for FHN model there is a grouping of kicks 
around the region where P2 is roughly zero. This means that kicks will not be felt as 
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Pi(6,0) 



-0.5 



-1 



Morris-Lecar 



0.2 



0.4 



0.6 . 0.8 



Fig. 10 The blue curves show the change in 9 under the action of a pulsatile kick in v, whilst the red 
dashed curves show the change in p under the same kick. The top plot is for the FHN model, whilst the 
bottom plot is for the ML model. We evaluate the effect of the kicks at p n = 0, where we observe the 
largest changes in 9 under the action of kicks 



much in the p variable, and so trajectories here do not get kicked far from cycle. This 
helps explain why it is more difficult to generate chaotic responses in the FHN model. 

After transients, we observe a 1 : 1 phase-locked state for the FHN model. For 
a phase-locked state, small perturbations will ultimately decay as the perturbed tra- 
jectories also end up at the phase-locked state after some transient behaviour. This 
results in a negative largest Lyapunov exponent of —0.0515. We note the sharply 
peaked distribution of kick phases, which is to be expected for discrete-time systems 
possessing a negative largest Lyapunov exponent, since such systems tend to have 
sinks in this case. The phase-locked state here occurs where P2 is small, suggesting 
that trajectories stay close to the limit cycle. Since kicks do not move trajectories 
away from cycle, there is no possibility of folding, and hence no chaotic behaviour. 
For the ML model, we observe chaotic dynamics around a strange attractor, where 
small perturbations can grow, leading to a positive largest Lyapunov exponent of 
0.6738. This time, the kicks are distributed fairly uniformly across 0, and so, some 
kicks will take trajectories away from the limit cycle, thus leading to shear-induced 
folding and chaotic behaviour. 



4 Discussion 

In this paper, we have used the notion of a moving orthonormal coordinate system 
around a limit cycle to study dynamics in a neighbourhood around it. This phase- 
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o.o 0.5 e n 1.0 o.o 0.5 e n 1.0 



Fig. 11 Panel (a) shows successive iterates of 6 in system (22)-(23) with functions P\2 taken from 
the FHN model, whilst panel (b) presents the same iterates but for functions P\2 from the ML model. 
Panel (c) shows P\ 2 for the FHN model, where the bold blue line is P\ and the red dashed line is P2 . 
Superimposed on this panel is a histogram displaying how kicks are distributed in terms of 6 alone. Panel 
(d) shows the same information, except this time for forcing functions from the ML model. Parameter 
values are a = 3.0, X = 0.1, s = 0.1, and T = 2.0 



amplitude coordinate system can be constructed for any given ODE system support- 
ing a limit cycle. A clear advantage of the transformed description over the original 
one is that it allows us to gain insight into the effect of time dependent perturba- 
tions, using the notion of shear, as we have illustrated by performing case studies of 
popular neural models, in two and higher dimensions. Whilst this coordinate trans- 
formation does not result in any reduction in dimensionality in the system, as is the 
case with classical phase reduction techniques, it opens up avenues for moving away 
from the weak coupling limit, where s —> 0. Importantly, it emphasises the role of the 
two functions P\(0, p) and P2(0) that provide more information about inputs to the 
system than the iPRC alone. It has been demonstrated that moderately small pertur- 
bations can exert remarkable influence on dynamics in the presence of other invariant 
structures [16], which cannot be captured by a phase only description. In addition, 
small perturbations can accumulate if the timescale of the perturbation is shorter than 
the timescale of attraction back to the limit cycle. This should be given particular 
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consideration in the analysis of neural systems, where oscillators may be connected 
to thousands of other units, so that small inputs can quickly accumulate. 

One natural extension of this work is to move beyond the theory of weakly cou- 
pled oscillators to develop a framework for describing neural systems as networks of 
phase-amplitude units. This has previously been considered for the case of weakly 
coupled weakly dissipative networks of non-linear planar oscillators (modelled by 
small dissipative perturbations of a Hamiltonian oscillator) [33-35]. It would be in- 
teresting to develop these ideas and obtain network descriptions of the following type: 

0i = 1 + fx (0i , Pi ) + J2 W iJ #1 (°i ' °J ' Pi ' Pj ) ' ( 26 > 

Pi = A(0i)pi +J2 w iJ H 2(0i>0j> Pi> Pj)> ( 27 ) 
j 

with an appropriate identification of the interaction functions H\ 2 in terms of the 
biological interaction between neurons and the single neuron functions P\ 2. Such 
phase-amplitude network models are ideally suited to describing the behaviour of the 
mean-field signal in networks of strongly gap junction coupled ML neurons [36, 37], 
which is known to vary because individual neurons make transitions between cycles 
of different amplitudes. Moreover, in the same network weakly coupled oscillator 
theory fails to explain how the synchronous state can stabilise with increasing cou- 
pling strength (predicting that it is always unstable), as observed numerically. All of 
the above are topics of ongoing research and will be reported upon elsewhere. 
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Appendix A: Derivation of the Transformed Dynamical System 

Starting from 

x = f(x)+eg(x,t), (28) 
we make the transformation x(t) = u(0(t)) + £ (0 (t)) p (t) , giving 

0 + S(0)P = f("(0) + mp) + eg(u(e + ?(0)p, t). (29) 

We proceed by projecting (29) onto §(#), using (1). The left-hand side of (29) now 
reads: 

du 
dO 



du(0) d?(0) ' 

p 

dO dO H 



+ q — P 
S d0 H 



dO 
df ' 



(30) 
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where § T denotes the transpose of § and the right-hand side of (29) becomes 



f f(u + i;p) + efg(u + i;p) 





du 






dO 


+ ? — P 
S d<9 _ 



+ $ T f(u + t;p) - ff{u) - ^ T ^- P + et T g(u + ; P ,t). 



(31) 



Thus, 



where 



e = i + M6, P ) + sh T (0, p) g (u{9) + op, 0. 



(32) 



and 



h(d,p) = 



Me, P ) = -h T (e, p)^p{0) + h T (e, P )[f(u + ; P ) - /(«)]. 

Upon projecting both sides of (29) onto £(#), the left-hand side reads 



(33) 



(34) 



du dt; 
dd + dd 



de dp 
dt dt 



T dt; de dp 

= c — P — — 

s d9 H dt dt 

= ? t ^p[i + Me, P ) + sh T (e, p) g (u + kp, 0] + %, 

de dt 



(35) 



whilst the right-hand side becomes 

f f{u + ;p) + s; J g(u + ;p,t) 
= -S T f(u) + ; r D/C - ; T Dft; + S T f(u + t;p) + sg(u + ;p,t), (36) 

since t, T f(u) = t} du/de = 0 and where Df denotes the Jacobian of /. Putting 
together the previous two equations yields 



p = A(e)p + Mo,p) + sr 

where 



t &m urn 
dO 



g(u(e) + t(e)p,t), (37) 



A(e) = s J 



de 



MO, p) = S T ^Ph +$ T [f(u + ?p) - /(«) - D/?p]. 



(38) 



(39) 
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It may be easily seen that f\(0, p) = O(p) as p -> 0 and that fi(0, 0) = 0 and 
9/2(0, 0)/3p = 0. Overall, combining (32) and (37) we arrive at the transformed sys- 
tem: 



! = 1 + MO, p) + sh T (0, p)g(u(0) + f (0, OP, 0' 



p = A(0)p + / 2 (0,p) + ^ i 



d£(0) 
In -^-ph(0,p) 



g(u(0) + W)p,t). 



(40) 



In order to evaluate the functions f\, fi, and A for models with dimension larger 
than two, we need to calculate d£ /d0 . Defining by Yi (0), the direction angles of § (0), 
we have that 



Ki = e i 



cos Yi 
1 + cos yi 
ei-m 



(«?i+§(0)), i=2,3,...,/i, 



(41) 



where the index / denotes the column entry of £ and x • y denotes the dot product 
between vectors x and y. Defining 



Ui (0) = 



and 



where j denotes the row index, we have 



dwi 



w j 



dui 
~d0 



dO dO 
By the quotient rule for vectors we find that 

dui (1 + ei)(et ■ (d£/d0)) - (e, • $(0))(ei(d£/d0)) 



d<9 



and that 



Overall, we have that 



d9 



dO d6 ' 



(42) 



(43) 



(44) 



(45) 



(46) 



(1 + eiKet ■ ~ (ei ■ £(0))(«i(d£/d0))\ 



(47) 
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Appendix B: Gallery of Models 

B.l Morris-Lecar 

The ML equations describe the interaction of membrane voltage with just two ionic 
currents: Ca 2+ and K + . Membrane ion channels are selective for specific types of 
ions; their dynamics are modelled here by the gating variable w and the auxiliary 
functions w^, r w , and m^. The latter have the form 

^ooO) = ^[1 +tanh(0 - vi)/v 2 )], r w (v) = l/cosh(0 - v 3 )/(2v 4 )), 

\ (48) 
Woc(v) = -[l +tanh((i; - u 3 )/u 4 )]. 

The function m^i;) models the action of fast voltage-gated calcium ion chan- 
nels; t>ca is the reversal (bias) potential for the calcium current and gca the cor- 
responding conductance. The functions z w (v) and u; 00(f) similarly describe the 
dynamics of slower-acting potassium channels, with its own reversal potential i>k 
and conductance #k- The constants ui ea k and gi ea k characterise the leakage current 
that is present even when the neuron is in a quiescent state. Parameter values are 
C = 20.0 uF/cm 2 , g\ = 2.0 mm ho/cm 2 , g K = 8.0 mmho/cm 2 , gca = 4.0 mmho/cm 2 , 
0 = 0.23, / = 39.5 uA/cm 2 , v\ = -60.0 mV, v K = -84.0 mV, u Ca = 120.0 mV, 
v\ = -1.2 mV, v 2 = 18.0 mV, v 3 = 12.0 mV, and v 4 = 17.4 mV. 



B.2 Reduced Connor-Stevens Model 



For the reduced CS model, we start with the full Hodgkin-Huxley model, with m, n, 
h as gating variables and use the method of equivalent potentials as treated in [26], 
giving rise to the following form for the function g: 



G(v, u) ■ 



dF 

~dh 



+ 



dF 

dn 



df dhoo(u) df driooiu) 



dhoo du 



dn c 



du 



(49) 



where dF/dh and dF /dn are evaluated 3th = hoo(u) and n = ftoo(w). For the gating 
variables (a,b), we have 



0oo 00 : 



r a (v) = 0.3632 



boo (v) = 



0.0761e^ +94 - 22 / 3L84 )\ 1/3 

1 _|_ e (u+1.17/28.93) 

1.158 



(50) 



1 +e O+55.96/20.12) ! 
4 



Tfc(v) = 1.24 + 



1 + e (u+53.3/14.54) 

2.678 



(51) 



1 4. e (u+50/16.027) ' 
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Parameter values are C = 1 uF/cm 2 , g\ = 0.3 mm ho/cm 2 , gK = 36.0 mmho/cm 2 , 
g a = 47.7 mmho/cm 2 , / = 35.0 u A/cm 2 , v 0 = 80.0 mV, v a = -75.0 mV, v K = 
-77.0 mV, vi = -54.4 mV, and u Na = 50.0 mV. 

B.3 FitzHugh-Nagumo Model 

The FHN model is a phenomenological model of spike generation, comprising of 2 
variables. The first represents the membrane potential and includes a cubic non- 
linearity, whilst the second variable is a gating variable, similar to w in the ML model, 
which may be thought of as a recovery variable. The system is 

l±v = v(a — v)(v — I) + I — w , w = v — bw, (52) 

where we use the following parameter values: fi = 0.05, a = 0.9, 7 = 1.1, and 
6 = 0.5. 
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